
# LIBRARIES ------

library(tidyverse)
library(readxl)
library(haven)


# DATA -------

# Result from models -----

ecuador_df <- read_rds("~/Downloads/ecuador_inpraise_df.rds")
india_kg_lpg <- read_rds("~/Downloads/india_kg_per_year_full_v3_oct16.rds")
kenya_kg_lpg <- read_rds("~/Downloads/kenya_kg_per_year_full_oct23.rds")


# Social cost of carbon ------
co2_scc <- read_rds("~/Downloads/1tco2_damages_1980_2020.rds") %>% 
  rename(Year = emission_year,
         scc = total_damages_dr2) %>% 
  add_row(Year = 1979,
          scc = 348) %>% 
  arrange(Year) %>% 
  mutate(
    scc_ratio = scc / 202.8917 # scc in 2020
  ) %>% 
  mutate(
    scc_low = 100*scc_ratio,
    scc_high = 1000*scc_ratio
  ) 

ggplot(co2_scc, aes(x=Year, y=scc, ymin=scc_low, ymax=scc_high)) + 
  geom_line() + 
  geom_line(data=co2_scc, aes(x=Year, y=scc_low), linetype="dotted") + 
  geom_line(data=co2_scc, aes(x=Year, y=scc_high), linetype="dotted") + 
  ggtitle("Social cost of carbon for 1 tCO2 emitted in a given year") + 
  ylab("USD") + coord_cartesian(ylim=c(0, 2000)) + theme_bw()

# ECUADOR -----

# countrywide: ~1970s = 4 megatonnes (million tonnes)
# countrywide: ~2019 = 41 megatonnes (million tonnes)

ecuador_co2e <- 
  ecuador_df %>% 
  left_join(
    co2_scc
  ) %>% 
  mutate(
    biomass_mj = polluting_pop * 7.23 * 365,
    biomass_mj_alt = polluting_pop_ctrfact_20yr * 7.23 * 365,
    
    lpg_mj = cf_pop * 7.23 * 365,
    lpg_mj_alt = cf_pop_ctrfact_20yr * 7.23 *365
  ) %>% 
  mutate(
    co2e_lpg = lpg_mj*173,
    co2e_lpg_ctrfact = lpg_mj_alt*173,
    
    co2e_biomass = biomass_mj*290,
    co2e_biomass_ctrfact =biomass_mj_alt*290,
  )  %>% 
  mutate(
    co2e = (co2e_lpg + co2e_biomass) / 1000000000000, # grams to megatonnes
    co2e_ctrfact = (co2e_lpg_ctrfact + co2e_biomass_ctrfact) / 1000000000000
  ) %>% 
  mutate(
    co2e_difference = co2e - co2e_ctrfact
  ) %>% 
  mutate(
    co2e_usd = co2e_difference * scc * 1e6,
    co2e_usd_low = co2e_difference * scc_low * 1e6,
    co2e_usd_high = co2e_difference * scc_high * 1e6,
  )

sum(ecuador_co2e$co2e_usd, na.rm=T)
sum(ecuador_co2e$co2e_difference, na.rm=T)
sum(ecuador_co2e$co2e, na.rm=T)
(abs(sum(ecuador_co2e$co2e_difference, na.rm=T)) - 
    sum(ecuador_co2e$co2e, na.rm=T)) / 
  sum(ecuador_co2e$co2e, na.rm=T)


# INDIA -----

# cost of subsidy 
india_subsidy_costs <- 
  india_kg_pm_scenario_1_yr_long %>%
  mutate(
    diff_900 = pred_1100_kg - pred_900_kg,
    diff_700 = pred_1100_kg - pred_700_kg,
    diff_550 = pred_1100_kg - pred_550_kg
  ) %>% 
  mutate(
    subidy_900_usd_total = (.93-.76) * diff_900,
    subidy_700_usd_total = (.93-.59) * diff_900,
    subidy_550_usd_total = (.93-.47) * diff_550
  ) %>% 
  mutate(
    subidy_900_usd_total_disc = subidy_900_usd_total / (1 + 0.09)^(Year-2023),
    subidy_700_usd_total_disc = subidy_700_usd_total / (1 + 0.09)^(Year-2023),
    subidy_550_usd_total_disc = subidy_550_usd_total / (1 + 0.09)^(Year-2023)
  ) %>% 
  group_by(model) %>% 
  mutate(
    subidy_900_usd_total_disc_cumsum = cumsum(subidy_900_usd_total_disc),
    subidy_700_usd_total_disc_cumsum = cumsum(subidy_700_usd_total_disc),
    subidy_550_usd_total_disc_cumsum = cumsum(subidy_550_usd_total_disc)
  ) %>% 
  group_by(Year) %>% 
  summarise( 
    yr_subidy_900_usd = mean(subidy_900_usd_total_disc, na.rm=T),
    yr_subidy_700_usd = mean(subidy_700_usd_total_disc, na.rm=T),
    yr_subidy_550_usd = mean(subidy_550_usd_total_disc, na.rm=T),
    
    subidy_900_usd = mean(subidy_900_usd_total_disc_cumsum, na.rm=T),
    subidy_700_usd = mean(subidy_700_usd_total_disc_cumsum, na.rm=T),
    subidy_550_usd = mean(subidy_550_usd_total_disc_cumsum, na.rm=T)
  ) %>% 
  pivot_longer(-Year) %>% 
  mutate(
    price = str_sub(name, 8, nchar(name))
  ) 


3.7*96000000*14.2 # mean refills * total HHs * 14.2kg per cyl = 5 043 840 000 off by factor of 10
 # 5 043 840 000 # back of envelope here

# 51 817 195 817 # from our models 

3.809e+10 / (3.7*96000000*14.2)

# India = 2.7 billion tonnes CO2e each year = 2700 megatonnes per year
# india_kg_lpg
india_co2e <- 
  india_kg_lpg %>% 
  mutate(
    lpg_observed_mj = kg_lpg_observed / 10 * 45, # 550 INR
    lpg_1_mj = kg_lpg_1100 / 10 * 45, # 1100 INR
    lpg_2_mj = kg_lpg_1 / 10 * 45, # variable INR
    # lpg_3_mj = kg_lpg_3 / 10 * 45
  ) %>% 
  mutate(
    lpg_observed_mj_deliv = lpg_observed_mj * .5,
    lpg_1_mj_deliv = lpg_1_mj * .5,
    lpg_2_mj_deliv = lpg_2_mj * .5,
    # lpg_3_mj_deliv = lpg_3_mj * .5
  ) %>% 
  mutate(
    delta_lpg_1_mj_deliv = lpg_1_mj_deliv - lpg_observed_mj_deliv,
    delta_lpg_2_mj_deliv = lpg_1_mj_deliv - lpg_2_mj_deliv,
    # delta_lpg_3_mj_deliv = lpg_3_mj_deliv - lpg_observed_mj_deliv
  ) %>% 
  mutate(
    delta_biomass_mj_1_deliv = delta_lpg_1_mj_deliv * (.5/.15),
    delta_biomass_mj_2_deliv = delta_lpg_2_mj_deliv * (.5/.15),
    # delta_biomass_mj_3_deliv = delta_lpg_3_mj_deliv * (.5/.15)
  ) %>% 
  mutate(
    delta_lpg_co2e_1 = delta_lpg_1_mj_deliv * 173,
    delta_biomass_co2e_1 = delta_biomass_mj_1_deliv * 275,
    
    delta_lpg_co2e_2 = delta_lpg_2_mj_deliv * 173,
    delta_biomass_co2e_2 = delta_biomass_mj_2_deliv * 275,
    
    # delta_lpg_co2e_3 = delta_lpg_3_mj_deliv * 173,
    # delta_biomass_co2e_3 = delta_biomass_mj_1_deliv * 275,
  )  %>% 
  mutate(
    delta_co2e1 = (delta_lpg_co2e_1 + delta_biomass_co2e_1)/ 1e12,
    delta_co2e2 = (delta_lpg_co2e_2 + delta_biomass_co2e_2)/ 1e12,
    # delta_co2e3 = (delta_lpg_co2e_3 + delta_biomass_co2e_3)/ 1e12
  )%>% 
  mutate(
    delta_co2e1_usd = delta_co2e1 * 203, # in millions
    delta_co2e2_usd = delta_co2e2 * 203,
    
    delta_co2e1_usd_low = delta_co2e1 * 100, # in millions
    delta_co2e2_usd_low = delta_co2e2 * 100,
    
    delta_co2e1_usd_high = delta_co2e1 * 1000, # in millions
    delta_co2e2_usd_high = delta_co2e2 * 1000,
    # delta_co2e3_usd = delta_co2e3 * 203
  ) %>% 
  mutate(
    delta_co2e1_usd_discounted = delta_co2e1_usd / (1 + 0.09)^(Year-2023),
    delta_co2e2_usd_discounted = delta_co2e2_usd / (1 + 0.09)^(Year-2023),
    
    delta_co2e1_usd_discounted_low  = delta_co2e1_usd_low  / (1 + 0.09)^(Year-2023),
    delta_co2e2_usd_discounted_low  = delta_co2e2_usd_low  / (1 + 0.09)^(Year-2023),
    
    delta_co2e1_usd_discounted_high = delta_co2e1_usd_high / (1 + 0.09)^(Year-2023),
    delta_co2e2_usd_discounted_high = delta_co2e2_usd_high / (1 + 0.09)^(Year-2023),
    # delta_co2e3_usd_discounted = delta_co2e3_usd / (1 + 0.05)^(Year-2023)
  ) 

india_co2e_summary <- 
  india_co2e %>% 
  group_by(
    model, price
  ) %>% 
  summarize(
    delta_co2e1 = sum(delta_co2e1, na.rm=T),
    delta_co2e2 = sum(delta_co2e2, na.rm=T),
    # delta_co2e3 = sum(delta_co2e3, na.rm=T),
    
    delta_co2e1_usd_discounted = sum(delta_co2e1_usd_discounted, na.rm=T),
    delta_co2e2_usd_discounted = sum(delta_co2e2_usd_discounted, na.rm=T),
    
    delta_co2e1_usd_discounted_low = sum(delta_co2e1_usd_discounted_low, na.rm=T),
    delta_co2e2_usd_discounted_low = sum(delta_co2e2_usd_discounted_low, na.rm=T),
    
    delta_co2e1_usd_discounted_high = sum(delta_co2e1_usd_discounted_high, na.rm=T),
    delta_co2e2_usd_discounted_high = sum(delta_co2e2_usd_discounted_high, na.rm=T),
    # delta_co2e3_usd_discounted = sum(delta_co2e3_usd_discounted, na.rm=T),
  ) %>% 
  group_by(
    price
  ) %>% 
  summarize(
    delta_co2e1_mean = mean(delta_co2e1, na.rm=T),
    delta_co2e2_mean = mean(delta_co2e2, na.rm=T),
    # delta_co2e3_mean = mean(delta_co2e3, na.rm=T),
    
    delta_co2e1_median = median(delta_co2e1, na.rm=T),
    delta_co2e2_median = median(delta_co2e2, na.rm=T),
    # delta_co2e3_median = median(delta_co2e3, na.rm=T),
    
    delta_co2e1_usd_discounted_mean = mean(delta_co2e1_usd_discounted, na.rm=T),
    delta_co2e2_usd_discounted_mean = mean(delta_co2e2_usd_discounted, na.rm=T),
    # delta_co2e3_usd_discounted_mean = mean(delta_co2e3_usd_discounted, na.rm=T),
    
    delta_co2e1_usd_discounted_median = median(delta_co2e1_usd_discounted, na.rm=T),
    delta_co2e2_usd_discounted_median = median(delta_co2e2_usd_discounted, na.rm=T),
    # delta_co2e3_usd_discounted_median = median(delta_co2e3_usd_discounted, na.rm=T)
    
    delta_co2e1_usd_discounted_low_mean = mean(delta_co2e1_usd_discounted_low, na.rm=T),
    delta_co2e2_usd_discounted_low_mean = mean(delta_co2e2_usd_discounted_low, na.rm=T),
    # delta_co2e3_usd_discounted_low_mean = mean(delta_co2e3_usd_discounted_low, na.rm=T),
    
    delta_co2e1_usd_discounted_low_median = median(delta_co2e1_usd_discounted_low, na.rm=T),
    delta_co2e2_usd_discounted_low_median = median(delta_co2e2_usd_discounted_low, na.rm=T),
    # delta_co2e3_usd_discounted_low_median = median(delta_co2e3_usd_discounted_low, na.rm=T)
    
    delta_co2e1_usd_discounted_high_mean = mean(delta_co2e1_usd_discounted_high, na.rm=T),
    delta_co2e2_usd_discounted_high_mean = mean(delta_co2e2_usd_discounted_high, na.rm=T),
    # delta_co2e3_usd_discounted_high_mean = mean(delta_co2e3_usd_discounted_high, na.rm=T),
    
    delta_co2e1_usd_discounted_high_median = median(delta_co2e1_usd_discounted_high, na.rm=T),
    delta_co2e2_usd_discounted_high_median = median(delta_co2e2_usd_discounted_high, na.rm=T),
    # delta_co2e3_usd_discounted_high_median = median(delta_co2e3_usd_discounted_high, na.rm=T)
  ) %>% 
  pivot_longer(-price)

india_co2e_summary

view(india_co2e_summary)




# india_co2e_alt <- 
#   india_kg_lpg %>% 
#   mutate(
#     kg_lpg_observed_mj = kg_lpg_observed * 45, # calorific value of 1kg LPG = 45 MJ / kg
#     kg_lpg_1_mj = kg_lpg_1 * 45,
#     kg_lpg_2_mj = kg_lpg_2 * 45,
#     kg_lpg_3_mj = kg_lpg_3 * 45
#   ) %>% 
#   mutate(
#     delta_lpg_mj_1 = kg_lpg_1_mj - kg_lpg_observed_mj,
#     delta_lpg_mj_2 = kg_lpg_2_mj - kg_lpg_observed_mj,
#     delta_lpg_mj_3 = kg_lpg_3_mj - kg_lpg_observed_mj
#   ) %>% 
#   mutate(
#     delta_lpg_kg_1 = kg_lpg_1 - kg_lpg_observed,
#     delta_lpg_kg_2 = kg_lpg_2 - kg_lpg_observed,
#     delta_lpg_kg_3 = kg_lpg_3 - kg_lpg_observed
#   ) %>% 
#   mutate(
#     biomass_mj_1 = delta_lpg_mj_1 * (.5/.15), # MJ Biomass = MJ LPG combusted * (50% efficiency of LPG stove / 15% efficiency of biomass stove)
#     biomass_mj_2 = delta_lpg_mj_2 * (.5/.15), # result = MJ delivered
#     biomass_mj_3 = delta_lpg_mj_3 * (.5/.15)
#   ) %>% 
#   mutate(
#     delta_lpg_co2e_1 = delta_lpg_kg_1*3725,
#     biomass_co2e_1 = biomass_mj_1* (.27) * 4200, # fnRB times CO2 only!
#     
#     delta_lpg_co2e_2 = delta_lpg_kg_2*3725,
#     biomass_co2e_2 = biomass_mj_2* (.27) * 4200,
#     
#     delta_lpg_co2e_3 = delta_lpg_kg_3*3725,
#     biomass_co2e_3 = biomass_mj_3* (.27) * 4200
#   )  %>% 
#   mutate(
#     delta_co2e1 = (delta_lpg_co2e_1 + biomass_co2e_1)/ 1000000000000,
#     delta_co2e2 = (delta_lpg_co2e_2 + biomass_co2e_2)/ 1000000000000,
#     delta_co2e3 = (delta_lpg_co2e_3 + biomass_co2e_3)/ 1000000000000
#   )%>% 
#   mutate(
#     delta_co2e1_usd = delta_co2e1 * 203,
#     delta_co2e2_usd = delta_co2e2 * 203,
#     delta_co2e3_usd = delta_co2e3 * 203
#   ) %>% 
#   mutate(
#     delta_co2e1_usd_discounted = delta_co2e1_usd / (1 + 0.05)^(Year-2023),
#     delta_co2e2_usd_discounted = delta_co2e2_usd / (1 + 0.05)^(Year-2023),
#     delta_co2e3_usd_discounted = delta_co2e3_usd / (1 + 0.05)^(Year-2023)
#   ) 

# india_co2e_summary <- 
#   india_co2e %>% 
#   group_by(
#     model, price
#   ) %>% 
#   summarize(
#     delta_co2e1 = sum(delta_co2e1, na.rm=T),
#     delta_co2e2 = sum(delta_co2e2, na.rm=T),
#     delta_co2e3 = sum(delta_co2e3, na.rm=T),
#     
#     delta_co2e1_usd_discounted = sum(delta_co2e1_usd_discounted, na.rm=T),
#     delta_co2e2_usd_discounted = sum(delta_co2e2_usd_discounted, na.rm=T),
#     delta_co2e3_usd_discounted = sum(delta_co2e3_usd_discounted, na.rm=T),
#   ) %>% 
#   group_by(
#     price
#   ) %>% 
#   summarize(
#     delta_co2e1_mean = mean(delta_co2e1, na.rm=T),
#     delta_co2e2_mean = mean(delta_co2e2, na.rm=T),
#     delta_co2e3_mean = mean(delta_co2e3, na.rm=T),
#     
#     delta_co2e1_median = median(delta_co2e1, na.rm=T),
#     delta_co2e2_median = median(delta_co2e2, na.rm=T),
#     delta_co2e3_median = median(delta_co2e3, na.rm=T),
#     
#     delta_co2e1_usd_discounted_mean = mean(delta_co2e1_usd_discounted, na.rm=T),
#     delta_co2e2_usd_discounted_mean = mean(delta_co2e2_usd_discounted, na.rm=T),
#     delta_co2e3_usd_discounted_mean = mean(delta_co2e3_usd_discounted, na.rm=T),
#     
#     delta_co2e1_usd_discounted_median = median(delta_co2e1_usd_discounted, na.rm=T),
#     delta_co2e2_usd_discounted_median = median(delta_co2e2_usd_discounted, na.rm=T),
#     delta_co2e3_usd_discounted_median = median(delta_co2e3_usd_discounted, na.rm=T)
#   ) %>% 
#   pivot_longer(-price)

# india_co2e_summary
# view(india_co2e_summary)


# KENYA ------


# VAT-generated $$$ 
sum(kenya_kg_lpg$kg_lpg_consumed_observed) * .2768

kenya_co2e <- 
  kenya_kg_lpg %>%  
  mutate(
    kg_lpg_consumed_observed_mj = kg_lpg_consumed_observed / 1.6 * 45,
    kg_lpg_consumed_1_mj = kg_lpg_consumed_1  / 1.6* 45,
    kg_lpg_consumed_2_mj = kg_lpg_consumed_2 / 1.6 * 45,
    kg_lpg_consumed_3_mj = kg_lpg_consumed_3  / 1.6* 45
    ) %>% 
  mutate(
    delta_lpg_mj_1 = kg_lpg_consumed_1_mj - kg_lpg_consumed_observed_mj,
    delta_lpg_mj_2 = kg_lpg_consumed_2_mj - kg_lpg_consumed_observed_mj,
    delta_lpg_mj_3 = kg_lpg_consumed_3_mj - kg_lpg_consumed_observed_mj
  ) %>% 
  mutate(
    biomass_mj_1 = delta_lpg_mj_1 * (.5/.15),
    biomass_mj_2 = delta_lpg_mj_2 * (.5/.15),
    biomass_mj_3 = delta_lpg_mj_3 * (.5/.15)
  ) %>% 
  mutate(
    delta_lpg_co2e_1 = delta_lpg_mj_1*173,
    biomass_co2e_1 = biomass_mj_1*395,
    
    delta_lpg_co2e_2 = delta_lpg_mj_2*173,
    biomass_co2e_2 = biomass_mj_2*395,
    
    delta_lpg_co2e_3 = delta_lpg_mj_3*173,
    biomass_co2e_3 = biomass_mj_3*395
  )  %>% 
  mutate(
    delta_co2e1 = (delta_lpg_co2e_1 + biomass_co2e_1)/ 1000000000000,
    delta_co2e2 = (delta_lpg_co2e_2 + biomass_co2e_2)/ 1000000000000,
    delta_co2e3 = (delta_lpg_co2e_3 + biomass_co2e_3)/ 1000000000000
  )%>% 
  mutate(
    delta_co2e1_usd = delta_co2e1 * 203 * 1e6,
    delta_co2e2_usd = delta_co2e2 * 203 * 1e6,
    delta_co2e3_usd = delta_co2e3 * 203 * 1e6,
    
    
    delta_co2e1_usd_low = delta_co2e1 * 100 * 1e6,
    delta_co2e2_usd_low = delta_co2e2 * 100 * 1e6,
    delta_co2e3_usd_low = delta_co2e3 * 100 * 1e6,
    
    
    
    delta_co2e1_usd_high = delta_co2e1 * 1000 * 1e6,
    delta_co2e2_usd_high = delta_co2e2 * 1000 * 1e6,
    delta_co2e3_usd_high = delta_co2e3 * 1000 * 1e6,
    
    
  ) %>% 
  mutate(
    delta_co2e1_usd_discounted = delta_co2e1_usd / (1 + 0.05)^(Year-2023),
    delta_co2e2_usd_discounted = delta_co2e2_usd / (1 + 0.05)^(Year-2023),
    delta_co2e3_usd_discounted = delta_co2e3_usd / (1 + 0.05)^(Year-2023),
    
    delta_co2e1_usd_discounted_low = delta_co2e1_usd_low / (1 + 0.05)^(Year-2023),
    delta_co2e2_usd_discounted_low = delta_co2e2_usd_low / (1 + 0.05)^(Year-2023),
    delta_co2e3_usd_discounted_low = delta_co2e3_usd_low / (1 + 0.05)^(Year-2023),
    
    delta_co2e1_usd_discounted_high = delta_co2e1_usd_high / (1 + 0.05)^(Year-2023),
    delta_co2e2_usd_discounted_high = delta_co2e2_usd_high / (1 + 0.05)^(Year-2023),
    delta_co2e3_usd_discounted_high = delta_co2e3_usd_high / (1 + 0.05)^(Year-2023),
    
    
  ) 
  
kenya_co2e_summary <- 
    kenya_co2e %>% 
  group_by(
    model, percent_price_change
  ) %>% 
  summarize(
    delta_co2e1 = sum(delta_co2e1, na.rm=T),
    delta_co2e2 = sum(delta_co2e2, na.rm=T),
    delta_co2e3 = sum(delta_co2e3, na.rm=T),
    
    delta_co2e1_usd_discounted = sum(delta_co2e1_usd_discounted, na.rm=T),
    delta_co2e2_usd_discounted = sum(delta_co2e2_usd_discounted, na.rm=T),
    delta_co2e3_usd_discounted = sum(delta_co2e3_usd_discounted, na.rm=T),
    
    delta_co2e1_usd_discounted_low = sum(delta_co2e1_usd_discounted_low, na.rm=T),
    delta_co2e2_usd_discounted_low = sum(delta_co2e2_usd_discounted_low, na.rm=T),
    delta_co2e3_usd_discounted_low = sum(delta_co2e3_usd_discounted_low, na.rm=T),
    
    delta_co2e1_usd_discounted_high = sum(delta_co2e1_usd_discounted_high, na.rm=T),
    delta_co2e2_usd_discounted_high = sum(delta_co2e2_usd_discounted_high, na.rm=T),
    delta_co2e3_usd_discounted_high = sum(delta_co2e3_usd_discounted_high, na.rm=T),
  ) %>% 
  group_by(
    percent_price_change
  ) %>% 
  summarize(
    delta_co2e1_mean = mean(delta_co2e1, na.rm=T),
    delta_co2e2_mean = mean(delta_co2e2, na.rm=T),
    delta_co2e3_mean = mean(delta_co2e3, na.rm=T),
    
    delta_co2e1_median = median(delta_co2e1, na.rm=T),
    delta_co2e2_median = median(delta_co2e2, na.rm=T),
    delta_co2e3_median = median(delta_co2e3, na.rm=T),
    
    delta_co2e1_usd_discounted_mean = mean(delta_co2e1_usd_discounted, na.rm=T),
    delta_co2e2_usd_discounted_mean = mean(delta_co2e2_usd_discounted, na.rm=T),
    delta_co2e3_usd_discounted_mean = mean(delta_co2e3_usd_discounted, na.rm=T),
    
    delta_co2e1_usd_discounted_median = median(delta_co2e1_usd_discounted, na.rm=T),
    delta_co2e2_usd_discounted_median = median(delta_co2e2_usd_discounted, na.rm=T),
    delta_co2e3_usd_discounted_median = median(delta_co2e3_usd_discounted, na.rm=T),
    
    delta_co2e1_usd_discounted_low_mean = mean(delta_co2e1_usd_discounted_low, na.rm=T),
    delta_co2e2_usd_discounted_low_mean = mean(delta_co2e2_usd_discounted_low, na.rm=T),
    delta_co2e3_usd_discounted_low_mean = mean(delta_co2e3_usd_discounted_low, na.rm=T),
    
    delta_co2e1_usd_discounted_low_median = median(delta_co2e1_usd_discounted_low, na.rm=T),
    delta_co2e2_usd_discounted_low_median = median(delta_co2e2_usd_discounted_low, na.rm=T),
    delta_co2e3_usd_discounted_low_median = median(delta_co2e3_usd_discounted_low, na.rm=T),
    
    
    delta_co2e1_usd_discounted_high_mean = mean(delta_co2e1_usd_discounted_high, na.rm=T),
    delta_co2e2_usd_discounted_high_mean = mean(delta_co2e2_usd_discounted_high, na.rm=T),
    delta_co2e3_usd_discounted_high_mean = mean(delta_co2e3_usd_discounted_high, na.rm=T),
    
    delta_co2e1_usd_discounted_high_median = median(delta_co2e1_usd_discounted_high, na.rm=T),
    delta_co2e2_usd_discounted_high_median = median(delta_co2e2_usd_discounted_high, na.rm=T),
    delta_co2e3_usd_discounted_high_median = median(delta_co2e3_usd_discounted_high, na.rm=T)
  ) %>% 
  pivot_longer(-percent_price_change)

kenya_co2e_summary
view(kenya_co2e_summary)


